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Abstract 

Renewable energy resources can indisputably minimize the threat of global warming and climate change. 
However, they are intermittent and need buffer storage to bridge the time-gap between production (off peak) 
and demand peaks. Based on geologic and geochemical reasons, the North German Basin has a very large 
capacity for compressed air/gas energy storage CAES in porous saltwater aquifers and salt cavities. Replacing pore 
reservoir brine with CAES causes changes in physical properties (elastic moduli, density and electrical properties) 
and justify applications of integrative geophysical methods for monitoring this energy storage. Here we apply 
techniques of the elastic full waveform inversion FWI, electric resistivity tomography ERT and gravity to map and 
quantify a gradually saturated gas plume injected in a thin deep saline aquifer within the North German Basin. 
For this subsurface model scenario we generated different synthetic data sets without and with adding random 
noise in order to robust the applied techniques for the real field applications. Datasets are inverted by posing 
different constraints on the initial model. Results reveal principally the capability of the applied integrative 
geophysical approach to resolve the CAES targets (plume, host reservoir, and cap rock). Constrained inversion 
models of elastic FWI and ERT are even able to recover well the gradual gas desaturation with depth. The spatial 
parameters accurately recovered from each technique are applied in the adequate petrophysical equations to 
yield precise quantifications of gas saturations. Resulting models of gas saturations independently determined 
from elastic FWI and ERT techniques are in accordance with each other and with the input (true) saturation 
model. Moreover, the gravity technique show high sensitivity to the mass deficit resulting from the gas storage 
and can resolve saturations and temporal saturation changes down to +3% after reducing any shallow fluctuation 
such as that of groundwater table. 

Keywords: Renewable energy; Compressed air/gas energy storage (CAES); Elastic full waveform inversion (FWI); 
Electric resistivity tomography (ERT); Gravity method; Petrophysical rock parameters 



Introduction 

One unprecedented challenge facing the human being is 
the energy resources, and its coupling with global climate 
changes and warming from greenhouse gases (GHG). 
Mitigation of anthropogenic GHG, including C0 2 emis- 
sions in the terrestrial atmosphere demands developments 
of viable alternative of renewable energy resources includ- 
ing hydroelectric, biomass, solar, wind, marine (wave/ 
tides) and geothermal sources. Most of these sources pro- 
duce energy only when suitable weather conditions are 
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prevailing and not when energy directly demanded. These 
sources are intermittent and need buffer storage to bridge 
the time-gap (disparity) between off-peak production and 
demand peaks. The underground geology offers an ad- 
equate option for short- and long-term energy storage 
such as compressed air or gas energy storage, CAES (e.g., 
Crotogino et al. 2001; Succar and Williams 2008). The 
North German Basin delivers favourable conditions (geo- 
logical, geochemical) for underground space utilization 
and has a huge capacity for CAES in porous brine aquifers 
and salt caverns (natural and artificial). Advantages of re- 
newable energy storage are (1) balancing power demand 
and fluctuating renewable energy production, (2) bridging 
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temporal mismatch between renewable energy production 
(off-peaks) and demand (peaks), i.e., storing off-peak en- 
ergy supply to use it during peak demand periods, and (3) 
offering large buffer capacity to meet any disruptions in 
energy supply. 

Figure 1 shows how a renewable electricity system could 
supply actual electricity demand during one week in 
Minnesota (Makhijani et al. 2012). The daily electricity de- 
mand is constantiy changing while the surplus renewable 
generation is put into storage. The basic approach of 
CAES is as follows: When electricity generation is greater 
than demand, the energy surplus is used to compress air 
in the geostorage. When generation is less than demand, 
pressurized air is withdrawn from storage and used to 
drive an electricity turbine. 

In a geological gas storage in saline formations, the 
gas replaces the pore brine causing strong changes in 
elastic moduli, density and electric resistivity. These 
physical contrasts justify the application of integrative 
geophysical techniques for monitoring this geostorage. 
These include techniques of elastic full waveform inver- 
sion (FWI), electric resistivity tomography (ERT), gravity 
and electromagnetic induction (EMI) of time- (TEM) 
and frequency-domain (FEM). EMI techniques (ground 
and air based) are applied usually for monitoring shallow 
targets, e.g., leakages in groundwater. 

Since some years ago Germany practices a turnaround 
in the energy policy (German "Energiewende") and is cur- 
rently leading in the production of solar and wind energy 
(IEA 2013). Wind energy are produced mainly at the 
coastal areas (on-shore and off-shore) of North Germany 
which is characterized by high wind speeds. We started 
2012 an interdisciplinary joint research project ANGUS+ 
dealing with impacts of using geologic subsurface as a 
thermal, electric or material storage in context with alter- 
native energy resources (Bauer et al. 2013). This includes 
dimensioning, risk analyses and impact predictions as a 
base for future space planning of the subsurface. Our 
main task is to develop a geophysical monitoring strategy 



using integrative approach of geophysical techniques 
(FWI, ERT, EMI and gravity) on almost realistic scenarios 
in the North German Basin. 

We show here results of numerical simulations of elas- 
tic FWI, ERT and gravity techniques in mapping CAES 
reservoirs with a continuous gradual desaturation with 
depth. These simulations are applied for synthetic data 
(before and after adding 3% random noise error) and 
inverted using constrains on the initial models. 

Gas storages in the North German Basin 

For a more realistic modelling scenario we selected a 
synthetic site in the North German Basin for this nu- 
merical study of CAES in geologic formations. The 
model block of this site (29*28*5.5 km 3 ) consists of a 
thick succession of 14 sedimentary layers ranging in age 
from Permian to Tertiary (Figure 2, Baldschuhn et al. 
2001; Hese 2012). The succession shows nearly horizon- 
tal layering with a gentle anticline fold in its southern 
part. It shows an unconformity, where formations of 
the Lower cretaceous, Lias and Rhaet disappear within 
the anticline crest at the depth interval of 0.7-1.0 km. 
The succession includes two thin brine reservoirs of a 
porous sandstone (5-30 m thickness), namely the Rhaet 
(1-1.5 km depth) and Quickborn (2-2.5 km depth) for- 
mations. Both are potential pore reservoirs for CAES 
and only the shallow Rhaet formation is considered in 
this study. The compressed gas is injected in the upper- 
most formation part within the northern anticline limb. 
Here this very light gas (e.g., air density of 1.29 kg/m 3 at 
STP "standard temperature, 0°C, and pressure, 101325 Pa") 
replaces the dense brine of >1100 kg/m 3 (with total dis- 
solved solids, TDS, of >100 g/1 at this depth) in the pores 
of the sandstone reservoir. Accordingly, the compressed 
gas saturation approaches its maximum value direcdy 
below the injection level and decreases gradually with 
depth within the reservoir of the dipping anticline flank. 
This downward gradual gas desaturation may correspond 
well with the almost realistic CO2 plume scenario injected 
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Figure 1 Daily supply and demand with storage of renewable energy, 11-17 July 2007, Minnesota (Makhijani et al. 2012). 



al Hagrey et al. SpringerPlus 2014, 3:267 
http://www.springerplus.eom/content/3/1/267 



Page 3 of 16 




Figure 2 Study site within the North German Basin, (a) 3D block below the Quaternary-Tertiary overburden, (b) 2D stratigraphy section 
perpendicular to the main anticline structure, and (c) 2D saturations of compressed gas (S gl and 5 gl ) injected in the saline Rhaet reservoir of the 
northern limb, (partly from Hese 2012). 



in a deep saline aquifer (Graupner et al. 2011). A combined 
simulation of multiphase flow, transport and geochemical 
reactions in a brine reservoir shows that the gas phase 
saturation decreases with increasing distance away from 
the injection point. The limited lateral extension of the 
thin Rhaet reservoir leads to dominating the downward 
gas propagation, where the porosity and permeability de- 
crease generally with depth. For this reservoir we applied 
two saturation distributions {S g i and Sgi) gradually decreas- 
ing downwards of the range 0.60-0.27 and 0.80-0.36, re- 
spectively (Figure 2c). The applied gradual desaturation 

with depth (z) follows this function: S g = Sgoe' 5 ^' 2 "^ , with 
Sgo = maximum S g at z 0 = 1 km, 1 km > z < 1.4 km. Both of 
S g i and Sg2 approach their maximum values near to the 
anticline crest (at 1 km depth) and gradually decrease with 
depth according to this function approaching their mini- 
mum values at 1.4 km depth. 

The seismic forward problem, FWI and model 
parameterization 

The propagation of seismic waves in an isotropic elastic 
medium can be described by a system of coupled first 
order partial differential equations (Landau and Lifschitz 
1986) 
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where d denotes the density, v L the particle velocity, o,y 
the stress tensor, Ey the strain tensor, A and \a. the Lame 



parameters, 8y the Kronecker Delta, ft, the source 
terms for body and surface forces, respectively. 

Equation (1) is a general expression for the conserva- 
tion of momentum in a continuum. It is independent of 
the medium state - such as gas, fluid or solid. To de- 
scribe the behavior of the material correctly a relation- 
ship between the forces (stresses ffy) acting on the 
medium and the resulting deformation (strain £y) is re- 
quired. For small forces/deformations and an isotropic 
medium this relationship is linear (generalized Hooke's 
law) depending only on the distibution of two material 
parameters A and {A (Lame parameters). Assuming that 
these parameters are time-independent the stress-strain 
relationship can be replaced by the stress-strain-rate 
equation (2). For a given isotropic elastic medium equa- 
tions (1-3) can be solved numerically and therefore syn- 
thetic seismograms for any acquisition geometry 
calculated. Based on the solution of the seismic forward 
problem a high-resolution imaging concept called full 
waveform inversion (FWI) has been developed in the 
1980s by Tarantola (1986). Since then the FWI is signifi- 
cantly improved and applied to a wide range of field ap- 
plications (Virieux and Operto 2009). 

The aim of (spatial) FWI is to minimize the data resid- 
uals Su = u mod - u obs between the modelled data u mod 
and the field data u obs to deduce high resolution models 
of the elastic material parameters in the underground. 
To solve this nonlinear optimization problem an appro- 
priate objective function E has to be defined. Similar to 
Asnaashari et al. (2013a), we use the following objective 
function 



\8u\ \ 2 + Ai 1 1 W m (m-m prior ) \\ 2 ) 



(4) 



Where the term \ \ 



\8u\\i denotes the L 2 -norm of the 
data misfit and \ \\ W m {m-m pr i or ) \\i a weighted L 2 - 
norm of the difference between the model parameters m 
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and prior model information m pr i or used for model 
regularization. The parameter Ai balances the contribu- 
tions of the data misfit and the model regularization 
term, while the spatial variable weighting factor W m de- 
fines which parts of the model are updated during the 
inversion process. The spatial weighting of the model 
updates is crucial for a successful inversion, because 
near surface inversion artefacts can introduce artificial 
data residuals not present in the real time-lapse data and 
therefore lead to an increase of the nonlinearity of the 
inverse problem. Like Asnaashari et al. (2013b) the mag- 
nitudes of the spatial weighting factors are based on 
elastic reverse-time migration (RTM) results to restrict 
model updates to the storage formation. The objective 
function equation (4) can be minimized by iteratively 
updating the model parameters m n (P-wave velocity 
(V p ), S-wave velocity (V s ), density (d)) at iteration step n, 
starting with an initial background model m 0 using the 
Newton method (Nocedal and Wright 2006): 

m n+ i = m n -T n {H-J l G m ) n (5) 

with the gradient of the objective function equation (4) 
with respect to the elastic model parameters 

Gm = \ P^^J + A i w ™ {m-m prior ) (6) 

and H m the second derivative of the objective function 
(Hessian). An explicit calculation of the Hessian in the 
time-domain is computational very expensive. Therefore, 
we use the quasi-Newton L-BFGS (Limited-memory 
Broyden-Fletcher-Goldfarb-Shanno) technique (Nocedal 
and Wright 2006; Brossier 2011), where the product of 
the inverse Hessian H m with the gradient G m is itera- 
tively approximated by finite-differences. 
The effective calculation of the time-domain gradient 

directions r^^ J with the adjoint method for different 

model parameterizations are described in Tarantola 
(1986), Mora (1987), Shipp and Singh (2002) and Kohn 
et al. (2012). The step length r n is estimated by a line- 
search satisfying the Wolfe-conditions (Nocedal and 
Wright 2006) to assure a fast and accurate convergence 
of the L-BFGS algorithm. While Asnaashari et al. 
(2013a) introduce another regularization term to assure 
model smoothness, we apply a weak wavenumber do- 
main filter to the estimated search directions at every it- 
eration step for the same purpose. 

In elastic time-lapse FWI, the data residuals are modi- 
fied according to 

Su= ({u m ° d (t 1 )-u mod (t 0 )}-{u obs (t 1 )-u° bs (t 0 )}), (7) 

which denote the difference between the modelled and 
the field data at time steps t 0 (baseline model) and t\ 



(Denli and Huang 2009). This redefinition of the data re- 
siduals leads to a much stronger focusing of the model 
updates at reservoir level. 

Based on the distribution of the air within the stor- 
age formation with a maximum gas saturation of 80% 
(Figure 2), an elastic model of the underground before 
and after the CAES injection is built. The elastic proper- 
ties of the rock matrix (P-wave velocity V Pim , S-wave vel- 
ocity V Sim , density d m and porosity <£») are linked with the 
physical parameters of the fluid and gas phases based on 
realistic matrix parameters to derive effective medium 
parameters. 

The effective seismic velocities V p and V s and the bulk 
density db of the 100% brine saturated aquifer before the 
CAES injection can be described by the following aver- 



aging equations (Gassmann 1951): 

d b = {\-<$>)d m + (t>d br (8) 

H=fid= (I"®) Pm (9) 
Kb = K d + {l-{K d /K m )) 2 /(0/K w + {1-<1>)/K m -K d /K 2 m ) (10) 

K d = (l-O) K m (11) 

K w =K br (12) 

V p = sqrt(K b /d b ) (13) 

V s = sqrt(fi b /d b ) (14) 



Where d denotes the density, K the bulk moduli and ^ 
the shear moduli. The subscripts b,m,d,br,g and w mean 
the bulk, rock matrix, dry rock, brine phase, gas phase 
and wetting phase, respectively. For the multiphase flow 
case after the CAES injection, equations (8) and (12) are 
adapted to: 

d b = (l-O) d m + O ((1-Sg)d br + Sgdg) (15) 
K w = SgKg + (l-Sg) K br (16) 

Based on this rock model Figure 3 shows changes of 
the material parameters due to the gas injection. Overall 
the changes of the effective medium parameters within 
the storage formation are quite small with maximum 
value changes of -240 m/s for V p , +86 m/s for V s and 
-168 kg/m 3 for d b - The small variations of the S-wave 
velocity are due to density variations only, while the 
shear modulus remains constant. 

Synthetic data of a reflection seismic survey along the 
transect is computed by solving the 2D isotropic elastic 
equations of motion equations (1-3) with a time- 
domain finite-difference (FD) technique on a Cartesian 
grid (Virieux 1986; Holberg 1987; Levander 1988). The 
reflection seismic acquisition geometry consists of 500 
vertical component geophones located at the surface. 
For the synthetic dataset 100 shots using a vertical impact 
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Figure 3 Input (true) changes in elastic parameters due to the gas injection in the Rhaet reservoir below the study site as obtained 
from the rock model, (a) AV P , (b) AV S and (c) density (Ad b ). 



source are recorded. The source signature is a 20 Hz 
Ricker wavelet. Model dimensions along the transect are 
6 km length and 1.8 km depth. Using an 8 th order spatial 
FD operator, the model can be discretized with 600 x 180 
grid points with a spatial grid point distance of 10 m. The 
time is discretized using At = 1 ms, thus for a recording 
time of t = 4 s, approximately 4000 time steps are com- 
puted. A free surface boundary condition is assumed on 
top of the model, while Convolutional Perfectly Matched 
Layers, CPMLs (Komatitsch and Martin 2007) are used at 
all other boundaries. The synthetic seismic sections are 
the input data for the elastic FWI. Figure 4 shows the 
common shot-gather for shot 50 of the baseline model 
and the time-lapse data between ti and to amplified by a 
factor 10. Notice the strong spurious multiple reflections 
due to the free boundary condition present in the 
baseline- and time-lapse data. 

To test the robustness of the elastic FWI approach 
for real field data applications we also investigate the 
influence of noise and added Gaussian noise within the 
frequency range of the source signal (0-40 Hz) with a 



signal-to-noise ratio S/N = 100 (1% noise) using the 
SUADDNOISE program of Seismic Unix (Cohen and 
Stockwell 2008). The resulting shot gathers are shown 
in Figure 4b,d for the baseline and time-lapse data, 
respectively. 

ERT modelling and parameterization 

At first we introduce briefly the approach for optimized 
electrode arrays in boreholes applied here. Like surface 
surveys, ERT data acquisition between two borehole 
electrode arrays can be conducted in the tripotential 
quadrupole configurations a (CPPC, C = current elec- 
trode, P = potential electrode), |3 (CCPP) and y (CPCP). 
For an N collinear multi-electrode array, a whole com- 
prehensive data set consists of [N(N- 1)(N- 2)(N- 3)/8] 
independent non-reciprocal quadrupole configurations 
(Noel and Xu 1991). The effective comprehensive data 
set results from excluding the redundant configurations 
of less stable inversions from the whole set, i.e., y config- 
urations and those of very large geometric factors (Loke 
et al. 2010). Resulting comprehensive data set is still 
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Figure 4 Common shot-gather for shot 50. (a, b) Baseline model and (c, d) time-lapse data residual between t, and t 0 amplified by a factor 
10. (a, c) Noise-free data and (b, d) noisy data with S/N = 100. 
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huge, e.g., a pair of 32 borehole electrodes yield >10 6 
data points. It can map subsurface targets with the high- 
est possible resolution but at very long acquisition times 
(i.e., poor temporal resolution) and at high costs. There- 
fore, an optimization approach is based on the model 
resolution matrix and searches for electrode configura- 
tions that maximize the resolution of survey results (e.g., 
Stummer et al. 2004; Wilkinson et al. 2006). An opti- 
mized borehole data sets of practical sizes (15,000 data 
points) of only 1.5% of the comprehensive data set but 
with almost the same spatial resolution is generated in 
this study (e.g., al Hagrey 2012a). Comparative applica- 
tions of diverse configurations (standard and non- 
standard) show the superiority of the optimized array re- 
sults (al Hagrey 2012b). Moreover, assessing our opti- 
mized array by the technique of region of investigation 
index (ROI) showed that its inverted tomograms are best 
constrained by the data coverage in comparison to that 
of the other configurations (Oldenburg and Li 1999). All 
these confirm the effectiveness of our optimization ap- 
proach applied here to generate a practical optimized 
data set of high resolution. 

The lstorage Rhaet formation consists of a highly re- 
sistive matrix (e.g., sandstone) and conductive pore brine 
saturant. The bulk resistivity (p) resulting from the gas 
displacing the brine is predicted using Archie's law 
(Archie 1942): 



where pb r is the brine resistivity, <D the porosity, S g the 
gas saturation and a,m and n are Archie constants. The 
separate phases (matrix, brine and gas) are assumed 
without any interaction. 

The electrical conductivity of the storage formation is 
caused mainly by the electrolytes of its pore brine. In the 
North German Basin, temperature and pressure, and par- 
ticularly the salinity or total dissolved solids (TDS) increase 
with depth. Increasing both TDS and temperature causes 
a dramatic decrease in the resistivity (e.g., Arps 1953; 
Schlumberger 1985). The TDS rise increases the number of 
ions carrying electrical currents. The temperature rise in- 
creases the salt solubility and decreases the brine viscosity 
which in turn enhances the ion mobility. The pressure in- 
crease with depth, on the other hand, causes a slight in- 
crease in the resistivity due to the closure of cracks that are 
often filled with conductive fluids. However, this effect de- 
creases with increasing depth and is negligible at pressure 
>0.3 GPa (e.g., Brace et al. 1965). In formations of the 
North German Basin, the average vertical gradient (with 
depth) of brine salinities, temperature and pressure ap- 
proach 100 mg/L/m, 0.03°C/m and 22.6 kPa/m, respectively 
(e.g., Magri et al. 2009). 



The target layers host borehole electrode arrays at 620 m 
offsets within a depth range of 0.9-1.8 km. Each array con- 
sists of 32 electrodes at 20 m spacing. The electric resistiv- 
ity of the 2D models are parameterized by the bulk rock 
resistivity values calculated from Archie equation (17). 
Here we applied the values of 0.2 for <t>, 0.08 Qm for pt, r 
(corresponding to TDS ~ 100 g/1 at 1 km depth of the 
North German Basin) and 1, 2 and 2 for constants a,m 
and n of, respectively, as typical values for the sandstone 
aquifer. Values of gas saturation are calculated by a poten- 
tial function simulating their gradual damping with depth. 

A 2.5D forward and inverse ERT modelling is carried 
out using modern codes (RES2DMOD, RES3DMOD x 
64 and RES2DINV x 64) based on algorithms by e.g., 
Loke et al. (2003). The forward modelling code is ap- 
plied to generate synthetic data sets between each adja- 
cent pair of borehole electrode arrays (including inhole 
and crosshole) using optimized electrode configurations. 
These synthetic data sets are generated after gas injec- 
tion in the brine reservoir of Rhaet formation. The data 
quality (0.6% average simulation error) is confirmed by 
results of tests on a homogeneous model with a constant 
p value. The technique robustness in the field is realized 
by adding a random error of commonly 3% to data sets 
in addition to their forward simulation error of 0.6%. 

In the ERT inversions, diverse setup constraints (mainly 
regularizations) are applied. These include the minimization 
methods of least squares {L 2 ) or robust blocky normalization 
(Li), and initial models of a constant homogeneous 
resistivity or an approximate inverse model (e.g., 
Claerbout and Muir 1973). Two synthetic data sets (gen- 
erated before and after adding 3% random noise) are 
inverted with incorporating mapping data of subsurface 
stratigraphy from prior (seismic) surveys (see next sec- 
tions). Each of these two data sets was inverted twice, 
once by incorporating layer interfaces and once by fixing 
resistivity regions, both are outside the reservoir layer. 

Gravity modelling and parameterization 

Rock densities depend on the mineral composition, por- 
osity and its content (fluid and gas), pressure p, 
temperature (T), deformations etc. The stratigraphy and 
average densities applied in the 3D gravity modelling of 
the study site are based on borehole measurements and 
data base of Geotectonic Atlas of northwestern Germany 
(Inselmann 1985; Baldschuhn et al. 2001; Hese 2012). 
The injected compressed gas in saline reservoirs dis- 
places pore brine and causes a drop of the bulk density 
which in turn causes a decrease in the gravity compo- 
nents and gradients. The bulk reservoir density (db) of 
partial gas saturations is given by equation (15), whereas 
dry air density is given by: 

d g =p/(RT) (18) 
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where R is the gas constant = 287.058 J/(kgK) for dry air. 
For calculating air and bulk densities (equations 18 and 
15), we considered almost realistic values for p (average 
27.786 GPa) and T (average 44°C) prevailing at 1- 
1.4 km depth range within the Rhaet reservoir of the 
North German Basin. This T value results from average 
local T at sea level (8°C) plus T fraction caused by the 
geothermal gradient down to 1.2 km depth (36°C). For 
instance, the average brine and air densities within this 
depth range approach 1123 and 303 kg/m , respectively, 
and the study sandstone reservoir will suffer from a bulk 
density drop up to -126 and -168 kg/m 3 corresponding to 
air saturations of S g i for Sgi, respectively. Figure 5 shows 
the 3D distributions of the gas saturation {S g i for Sgi) and 
bulk density (resulting from the parameterization) in 
the brine reservoir of the anticline flank within the North 
German Basin. 

We used here the software IGMAS+ (Interactive 
Gravity and Magnetic Application System) designed for 
3D gravity, gravity gradient and magnetic modelling 
(e.g., Gotze and Lahmeyer 1988; Schmidt et al. 2011). 
The model is extrapolated outside the volume of interest 
in all directions (about two times the model length) to 
avoid any edge effects. We calculated the gravity field 
components (g z , g y and g z ) and gradients (g zx , g zy and g zz ) 
before and after CAES injection, respectively, as well as 
their difference (residual) anomalies (Ag^, Ag y and Ag z ). 
Here we show vertical component Ag z maps only which 
reflect the strongest anomalies with respect to CAES 
reservoirs. We will discuss these gravity anomalies 
resulting from the two saturations (S g i and Sgi) and the 
lower sensitivity boundary determined for the technique 
at all saturations, i.e., the least measurable gravity anom- 
aly determined by the modern micro-gravimeter accur- 
acy (3-5 uGal). 




Figure 5 3D distribution of gas saturation (S g ) and bulk density 
(d) in the brine reservoir of the anticline flank (see Figure 2c). 

Reference d of reservoir (i.e., at S g = 0) is 2154 kg/m 3 . 



Results of elastic time-lapse FWI 

The initial model for the time-lapse waveform inversion at 
each time-step is the true elastic medium model before 
the CAES injection. While this seems to be an overopti- 
mistic assumption, we want to focus this part of the study 
on the question if the elastic FWI is capable to reconstruct 
structures at the resolution limit at all. Later we will also 
investigate the impact of different errors in the baseline 
model on the elastic time-lapse FWI results. To reduce 
the nonlinearity of the multiparameter inversion problem 
a sequential frequency FWI approach is applied. There- 
fore, Butterworth-lowpass filters are applied to the source 
wavelet and field data with corner frequencies of 20 and 
40 Hz, respectively. To reduce the influence of multiple 
reflections, exponential time damping (Brossier et al. 
2009) of the time-lapse data after the first arrivals are ap- 
plied. The first arrivals are automatically picked with a 
STA/LTA picker for the initial model. The detailed elastic 
FWI workflow is described in Table 1. All material param- 
eters are simultaneously updated. 

The inversion results of the noise-free and noisy seismic 
time-lapse data (Figure 6) can be compared with the true 
changes of V p , V s and bulk density db in Figure 3. For the 
noise-free data (Figure 6, top panel) the shape of the 
CAES plume, the gradient in the seismic velocities and 
density could be reconstructed well. The magnitudes of 
the different material parameters also agree with the true 
model. The resolution limit of the elastic FWI is roughly 
one quarter to half of the minimum seismic wavelength. 
Using a maximum frequency of 60 Hz with a S-wave vel- 
ocity of 2074 m/s at reservoir level leads to a minimum 
seismic wavelength of 34 m and therefore a resolution 
limit of 9-17 m. This coincides with the FD model 
discretization errors of 10 m at the layer interfaces visible 
in the FWI results. Using the noisy time-lapse data 
(Figure 6, bottom panel) the resolution of AV P and Ad 
models is also quite good, beside some minor artefacts. 
Only the AV S model is affected by the Gaussian noise. Fur- 
ther tests (not shown here) with a S/N = 50 (2% noise) 
lead to even more dominant noise artefacts, but the gas 
plume is still visible. At S/N = 25 (4% noise) the gas plume 
vanishes within the noise artefacts. 

Gas quantification by elastic time-lapse FWI 

In most cases gas saturations are estimated from the 
P-wave velocity of elastic FWI results (Queisser and 

Table 1 FWI workflow: corner frequencies (/ c ) of 
Butterworth-lowpass filter for the sequential frequency 
inversion and time-damping coefficients 

f c [Hz] Time-damping coefficients [1/s] 

20 100,50,5,1 
40 1 00, 50, 5, 1 
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Figure 6 FWI results showing changes in seismic parameters, (a, b) P-wave velocity Al/ p , (c, d) S-wave velocity Al/ S and (e, f) density Ad for 
the noise-free and noisy dataset (S/N = 100), respectively (cf., Figure 3). 



Singh 2012) by inverting the Gassmann equation (10). 
This involves the estimation of a lot more or less un- 
known petrophysical parameters. Because the elastic 
multiparameter FWI delivers a density model, the inver- 
sion of equation (15) seems to be more appropriate. By 
assuming that the brine and gas densities remain con- 
stant, the following equation can be derived from equa- 
tion (15) for the CAES saturation: 



d t -do 

0{d g -dbr) 



(19) 



depending only on the change of the bulk density between 
time t (d t ) and the baseline model (do), the porosity, as 



well as brine and gas density. The correct values for por- 
osity and bulk density of the baseline model in the storage 
formation are used. In Figure 7, the estimated gas satur- 
ation changes from the elastic FWI density model (b, c) is 
compared with the true saturation (a) for the noise-free 
(b, d) and noisy data (c, e). Using the noise-free data, the 
extension of the CAES phase and the CAES saturation 
changes are well recovered. Due to the finite frequency 
content of the source wavelet, the sharp boundaries of the 
CAES plume cannot be resolved. Therefore, fictitious 
large absolute gas saturation errors of about ±10-15%, lo- 
cally up to 30%, occur at the boundaries (Figure 7d, e). 
The average error within the reservoir approaches about 
5%. With the density model from the noisy FWI result 
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Figure 7 Gas saturation after the CAES injection, (a) True saturations S g " ue , (b, c) S g estimated from the density model of the FWI result, 
and (d, e) absolute error for the noise-free data and the noisy data (S/N = 100), respectively. 
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absolute saturation errors increase to 30-60% at the 
plume boundaries and 5-30% within the plume. Noise ar- 
tefacts outside the plume can lead to errors of up to 20%. 
For larger S/N-ratios estimated gas saturations become 
unrealistically larger than 100%. 

Impact of errors in the baseline model on the elastic 
time-lapse FWI results 

So far we only investigated the resolution of elastic time- 
lapse FWI when the baseline model is perfecdy known. 
Figure 8 shows a workflow for the estimation of different 
baseline models with different kinds of errors. In the first 
step smooth versions of the true baseline seismic velocity 
models are generated (Figure 8a), while the density model 
is assumed to be constant. 



With the macro velocity models a pre-stack RTM can 
be applied to the seismic data to reconstruct the posi- 
tions of the layer interfaces (Figure 8b). In the next step 
the interfaces are manually picked (Figure 8c) in the 
zero-offset section and information about the V p , V s and 
db model at a borehole location are used to fill the layers 
with material parameters (Figure 8d, baseline model A). 
A correct localization of the layer interfaces can be 
quite complicated, especially in the gas storage layer. 
As a result velocity and density errors of up to 1000 m/s 
and 300 kg/m 3 , respectively, are introduced in the base- 
line model, even though the medium properties within 
the layers are lateral homogenous. Beside erroneous 
multiple reflections within the storage formation due to 
wrong interface positions, a different diffraction pattern 
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due to changes of model discretization errors on the 
Cartesian FD-grid become a problem. The elastic 
time-lapse FWI result for baseline model A is com- 
pared in Figure 9d-f with the result for the perfect 
baseline model (Figure 9a-c). While the position and 
shape of the gas plume can be reconstructed at least 
to some extent, the distribution of the material proper- 
ties within the gas plume are systematically overesti- 
mated. To improve the quality of baseline model A, it 
is smoothed (Figure 8e) with a subsequent application 
of spatial FWI to the baseline data (Figure 8f, baseline 
model B). This additional spatial FWI step reduces 
the errors in the velocity models at the layer inter- 
faces up to a factor two, but the errors are still quite 
large with ±500 m/s for V p , ±400 m/s for V s , and 
± 200 kg/m 3 for db- Additionally the smoothing introduces 
velocity and density errors within the layers. As a result 
the time-lapse FWI results of baseline model B are only 
slightly improved (Figure 9g-i) compared to baseline 
model A (Figure 9d-f). 

The evolution of the objective function E during the 
elastic time-lapse FWI is shown in Figure 10. Transitions 
between the different inversion stages are defined in 
Table 1. For the perfect baseline model without noise 
the objective function significantly decreases during the 
FWI, while none of such significant decrease is notable 
when data with a S/N = 100 is inverted. The number of 
iterations is also reduced. A comparable behaviour oc- 
curs for baseline model A and B. 



ERT inversion results 

As mentioned before the two synthetic data sets (generated 
before and after adding 3% random noise) are inverted with 
incorporating mapping data of subsurface stratigraphy from 
prior surveys (e.g. seismic, borehole logs). These two appar- 
ent resistivity data sets were inverted by incorporating layer 
interfaces and fixing resistivity regions, respectively, both 
are outside the reservoir layer. ERT data inversions recon- 
struct directly the true subsurface resistivity tomograms 
including the study gas plume of downward gradual desat- 
uration. This implies no model differencing, unlike seismic 
results showing the model difference before and after the 
gas injection. Of all differently independent inversions, 
every best-fitting tomogram shows least root mean square 
(rms)-errors of <0.5% and iteration number almost of 5, 
and is optimized with the L\ norm for sharp interfaces. 
Also this low rms-error value is explained by the good con- 
vergence of the synthetic data sets toward the final solution. 
This Li norm yields significantly more accurate results than 
L 2 norm, where the actual subsurface resistivity changes 
abruptly at sharp target boundaries (cf. Loke et al. 2013). It 
is more likely to suppose that considerable subsurface infor- 
mation is already available during monitoring from the de- 
tailed baseline survey and any other subsequent survey. 
The a priori incorporation of this mapping data in ERT in- 
versions minimizes the ambiguity of the solution and en- 
hances the resolution of results. Here each data set 
inversion is constrained by incorporating resistivity regions 
and interfaces outside the reservoir, respectively. 
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Figure 9 Influence of different baseline models on the elastic time-lapse FWI results, (a-c) Changes of Vp, Vs and density for the perfect 
baseline model, (d-f) baseline model A estimated from RTM results and borehole logs, and (g-i) baseline model B estimated from RTM results, 
borehole logs and spatial FWI. 
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Figure 10 Development of the objective function E (equation 1) during the time-lapse FWI for the perfect baseline/noise free case, 
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Moreover, we calculated the model resistivity differ- 
ence (residual) relative to the input (true) model to 
quantitatively evaluate the reliability of the recon- 
structed tomograms. This residual anomaly is a meas- 
ure for the reliability of ERT inversions. It expresses 
quantitatively the resolution of the applied technique 
for both of the spatial mapping capability and the re- 
covering resistivity amplitude. An accurate resistivity 
amplitude is essential for precise quantifying the injected 
gas phase saturation using Archie equation (17), see next 



section. This residual (Ap) between the corresponding 
pixels of the input or true {p trU e) and output (pert) 2D 
model is calculated by: 



Ptme Pert i 



(20) 



Figure 11 shows only the best-fitting tomograms re- 
constructed together with their corresponding residuals 
for the four cases of applying synthetic data before and 
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Figure 11 ERT modelling results, (a) Input or true (p true ) models, (b, c, f, g) output resistivity (fl ERT ) tomograms with rms-error, and (d, e, h, i) 

corresponding residual Ap models with their average values. Inversion is conducted using constrains of incorporating p region (b-e) and boundaries 
(f-i). These inversions are conducted for synthetic borehole ERT data (contaminated with 0.6% simulation errors) before and after adding 3% random 
noise error. The low residuals (average <9%) reflect well the advantage of constraining of ERT inversions by the seismic mapping results for an accurate 
reconstruction of the gas phase of downward gradual desaturation. The vertical black dots mark the borehole electrodes and the continuous lines 
mark the interfaces. 
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after adding 3% random noise errors, and constrained 
inversions of incorporating p regions and boundaries 
outside the reservoir layer. 

Obviously, the ERT inversion results of reconstructed 
tomograms visualize well the storage targets, particu- 
larly the gas plume of a downward gradual desaturation 
within the brine reservoir for all studied four cases 
(Figure 11). It is clear that constrained inversion models 
incorporating resistivity regions are better resolved than 
these incorporating boundaries. Also the addition of 3% 
random error to the synthetic data sets increases the 
misfit of rms-error values (between input and output 
response) by a factor of 5-8 but slightly decreases the 
mapping resolution as reflected by low residual Ap rise 
from 4.5 and 7.8% to 4.9 and 7.9%, respectively. Loke 
et al. (2013) found that the model inverted with the 
ii-norm is less sensitive to random noise compared 
with the L 2 -norm. 

Resulting residual Ap r tomograms as a measure for 
resolution confirm generally the good mapping capabil- 
ity of the applied constrained ERT technique, where all 
average Ap r have low values of 4.5-7.9% with minor de- 
viations (Figure 11). This Ap r distribution shows that the 
resulting model reliability (inverse of Ap r ) is least for the 
noisy data set inverted by incorporating boundaries and 
best for the data set without adding random noise 
inverted with fixing p regions. An accurate investigation 
of Ap r tomograms show that the resolution suppression 
due to the addition of 3% random errors to the gener- 
ated data (average Ap r increases from 4.5 to 4.9 for fix- 
ing p region, and from 7.7 to 7.9 for incorporating 
boundaries) is lower than that resulting from incorporat- 
ing boundaries instead of fixing p regions (average de- 
clines from 4.5 to 7.7 for data without random noise, 
and from 4.9 to 7.9% for noisy data). 

Constrained inversion tomograms (using any available 
subsurface data as an a priori information in the ERT in- 
version) show better resolution (inverse of Ap„ Figure 11) 
than their corresponding unconstrained or even partly 
constrained inversion tomograms (not shown here, see 
al Hagrey et al. 2013). However, these residual Ap maps 
still reflect the common smearing effects and artifacts of 
varying degrees of the ERT technique. This negative effect 
is particularly visible within the thin storage formation 
(Ap r up to ±20%) with spatially varying gas saturation (i.e., 
resistivity amplitude). On the other hand the non-varying 
(homogeneous) formations above (cap rock) and below 
(aquitard) the storage formation show a Ap r of almost ±0. 
Obviously almost all inversion uncertainties are related to 
the monitored gas phase within the host reservoir. Thus 
they deliver an error estimate of the injected gas quantities 
monitored by this constrained inversion technique. 

In conclusion, the ERT technique with permanently in- 
stalled borehole electrodes aims at mapping, monitoring 



and quantifying the gas volume injected into the saline 
aquifer at any time. Obviously, the resulting tomograms 
(Figure 11) fulfil well the spatial mapping and monitoring 
purposes. This permanent electrode installation helps to 
maximize the reliability of monitoring data. Modelled 
anomalies minimizes the background effect and thus max- 
imizes the time-varying response caused here by the 
injected gas quantity. The residuals (Figure 11) assess and 
prove the high reliability of the results including the quan- 
tification capability for the resistivity amplitude. These 
highly reliable resistivity amplitudes motivated us to derive 
the gas saturations (see next section). 

Gas quantification by ERT 

The gas phase saturation S g or S g ERT in a partially satu- 
rated reservoir medium is driven indirectly by applying 
the resistivity amplitudes recovered from ERT models 
(Figure 11) in the Archie equation rearranged as: 

Here, we assume no interaction between the three res- 
ervoir phases (solid matrix, brine and gas). Using this 
equation, saturation models of gas phase (S g ERT ) within 
the storage reservoir are calculated and plotted together 

with the corresponding input [true, Sg Ue ^J models and 

their absolute difference Us g = S^ ue -Sf T ^j in Figure 12 

for the four applied cases described before. It is clear 
that the estimated saturation here contains the uncer- 
tainty of the other parameters of a,m,n,<I> and pi, r in- 
cluded in this equation. Unlike this, the equation of the 
resistivity index (the resistivity ratio after and before the 
injection, e.g., Nakatsuka et al. 2010) eliminates parame- 
ters a,m and <t>, whereas parameters pb,- and n are as- 
sumed to stay unchanged with time. However, this index 
equation includes the uncertainties of both pre- and 
post-injection p models unlike our applied equation (21) 
which includes the uncertainty of the post injection 
model only. In the porous sandstone of the Rhaet reser- 
voir, we considered the common values of 1, 2 and 2 for 
the constants a,m and n, respectively (Archie 1942), 0.2 
for O, and 0.08 flm for pb r (as mentioned before). 

An investigation of resulting models shows that the 
gradual gas desaturations {S g ) with depth are well 
mapped within the thin storage formation. The results in 
Figure 12 show a satisfactory similarity between input 
(true) and corresponding output (reconstructed) models. 
The average difference AS g relative to s g '" put (Figure 12) 
approaches 15-21% for all studied four cases of Figure 11. 
This small difference confirms again the satisfactory reli- 
ability of the results. The gas distribution within the Rhaet 
reservoir formation shows a good similarity particularly in 
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Figure 12 Results of gas phase saturations in the storage formation derived from ERT models (of Figure 11). (a) True simulation (S g '" put ), 
(b, c, f, g) saturations (S g ERT ), and (d, e, h, i) corresponding differences (AS g ) with their average values. The gas phase is recovered satisfactorily 
everywhere as reflected by AS g values (<21%). The continuous lines mark reservoir boundaries. 



the region 1-1.4 km (corresponding to the real distri- 
bution) and is insignificant outside this region where its 
amplitude is below the average {ASg). Most recovered 
S g ERT values are lower than their corresponding input 
ones. This may be related to the smearing effects of the 
technique. This smearing influences negatively the 
modelled resistivity amplitude and reduces the resistiv- 
ity high of the gas plume sandwiched between the two 
resistivity lows of caprock and aquitard, respectively. 
High ASg values are concentrated mainly at reservoir 
interfaces and may be related to discretization errors. 
At the Cranfield site, Mississipi, Carrigan et al. (2013) 
found that CO2 saturations measured in monitoring 
wells are higher than the ERT-derived saturations al- 
though both show good spatial correlations. They 
added that ERT provides an integrated response from 
large volume, whereas gas sensors (dm penetration) 
provide point measurements and are sensitive to condi- 
tions near the well. 

3D gravity modelling results 

Results of the 3D gravity modelling technique using 
IGMAS+ program show a high sensitivity to the applied 
plume scenarios of gas phase injected in the pore Rhaet 
reservoir of the North German Basin (Figure 13). The 
negative anomaly amplitude of the vertical gravity compo- 
nent after subtraction the background before gas injec- 
tions (Ag z ) increases with increasing the gas saturation 
S g causing the mass deficit. The two saturations of S g i 



(60-27%) and S g2 (80-36%), gradually decreasing with 
depth, show negative anomalies Ag zl and Ag z2 , respect- 
ively, of similar shape, i.e., visually both are hardly dis- 
tinguishable (Figure 13a, b). The absolute Ag zl 
amplitude for S g i (134 uGal) is less than that for S g 2 for 
Sg2 (178 uGal). The difference in Ag z amplitude (i.e., 
double difference) between both saturations is visual- 
ized well in figure (Figure 13c) and is far higher (>10 
times) than the measurable accuracy of modern micro- 
gravimeters (±3-5 uGal). Fewer saturation changes are 
verified systematically at the whole range of saturations. 
We found that 1% change in saturations AS g yield 2.2 
uGal change in Ag z for our study gas plume. This im- 
plies that the technique can monitor a saturation 
change AS g down to 2.5-3% for the whole saturation 
range. This AS^ range results in a Ag z changes above 5 
uGal. Obviously, time-lapse data do not require many 
corrections (e.g., free-air, Bouguer, terrain) but temporal 
shallow changes (e.g. fluctuations of the water table) 
highly affect the gravity readings. Our gravity modeling of 
a water table at 10 m depth below the study site yields a 
measurable 5 uGal anomaly (micro-gravimeter accuracy) 
already for 0.5 m fluctuations only. Therefore, such fluctu- 
ations should be observed in wells to remove their effects 
from gravity readings. In conclusion, the 3D gravity mod- 
elling technique applied here is able to monitor time- 
lapses of gas saturations of S gl and Sg2 as well as any 
saturation changes (AS g ) down to ±3% within the whole 
range of saturation. 
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Figure 13 3D forward gravity modelling of gas storages in the anticline limb of Rhaet sandstone formation below the study site. 

(a) Vertical gravity anomalies Ag z (relative to background) for the gas reservoir with saturations 5 g] , (b) Ag z for 5 g2 , (c) their anomaly difference, 
and (d) least anomaly difference measurable by micro-gravimeters of >5 uGal resulting from a saturation change AS g of ±3%. Only the measurable 
anomaly part is visualized in c-d by blanking colour bar below 5 uGal (immeasurable limit). 



Discussion and conclusion 

Mitigation of anthropogenic green house gases demand 
developments of renewable energy resources. However, 
most renewable energy sources are intermittent and 
therefore need buffer storage to match public power 
supply and demand. In geologic storages, replacing pore 
brine with compressed gas energy within the reservoir 
formation causes changes in elastic and electric proper- 
ties and justify applications of integrative geophysical 
methods for monitoring this energy storage. We apply 
here elastic FWI, ERT and gravity techniques to map 
and quantify a thin gas plume of gradual downward de- 
saturation injected in a deep brine aquifer. These aimed 
tasks are real challenges for any singly applied geophys- 
ical monitoring technique. 

In this numerical study we simulate a nearly realistic 
storage scenario by considering: (1) the study site is 
chosen inside the North German Basin of favourable 
conditions for energy storages, (2) a storage model sce- 
nario is parameterized by real (published) data for this 
basin, (3) the gas phase plume is simulated with down- 
ward gradual desaturation similar to realistic cases, and 
(4) a common random noise level is added to the syn- 
thetic data to robust the technique for real field applica- 
tions. The aimed resolution is enhanced by applying: 
(1) an integrative approach of geophysical methods, (2) 
an optimized approach for data acquisitions with a data 
coverage constraining well the inversion model and maxi- 
mizing the resolution, and (3) constrained inversions to 
minimize interpretation ambiguities (by a priori use of 
available data, e.g. seismic, logs). 

Unlike classical travel-time based tomographic ap- 
proaches, the elastic FWI is capable to map the exten- 
sion of the thin gas plume of downward gradual 
desaturation using only reflection seismic data, if a very 
accurate background model for the seismic velocities 
and density can be estimated before the gas injection. 



Additionally the elastic FWI recovers the changes of iso- 
tropic elastic material parameters and density due to the 
gas injection and subsequent partial drainage of the 
aquifer. By using an appropriate rock model, changes of 
the gas saturation can be deduced from the elastic FWI 
results with an accuracy of 5-30% within the aquifer, de- 
pending on the amount of noise (S/N-ratio > 100) 
present in the data. Due to the finite frequency content 
of the source signal larger saturation errors up to 60% 
can occur at the boundaries of the gas plume. Density 
inversion artefacts outside the aquifer due to noise can 
lead to fictitious estimates of saturation variations with 
local maxima of 20%. For a S/N-ratio of 50 the shape of 
the gas plume is still visible, but estimations of the gas 
saturation become highly erroneous and a S/N-ratio of 
25 seem to be the detection limit for the gas plume. Er- 
rors in the elastic baseline model has a very substantial 
impact on the quality of the elastic time-lapse FWI re- 
sults. Picking errors in layer interfaces and inaccurate 
material parameters within the layers lead to results 
with an approximately correct shape and position of 
the gas plume, but overestimated wrong elastic material 
parameters within the gas plume. Therefore, the estima- 
tion of accurate elastic baseline models for a successful 
elastic multiparameter FWI, and a subsequent calculation 
of the gas saturation distribution, is the greatest challenge 
for real field applications. Smooth macro-velocity models 
based on Common-Reflection Surface (CRS) stacking 
(Mann 2002) and Normal-Incidence Points (NIP) -wave 
tomography (Duveneck 2004) for P- and SH-wave data 
combined with the intensive use of prior information from 
borehole logs seem to be the most promising approach. 

The applied constrained ERT inversion technique (tak- 
ing use from previous seismic mapping and well logs) is 
also able to accurately map storage targets (caprock, res- 
ervoir with thin gas and aquitard) in all four applied 
cases (resulting from inverting noise-free and noisy data 
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by incorporating layer boundaries and resistivity regions, 
respectively). The technique can even recover the gas 
plume of downward gradual desaturation with a good 
resolution. Also inversion models constrained by incorp- 
orating resistivity regions are better resolved than these 
constrained by incorporating boundaries, both are ap- 
plied outside the reservoir layer. The thin resistive gas 
plume is sandwiched between the two conductive layers 
of the overlying caprock and the underlying aquitard. 
Based on the equivalence principle, the resistance (p*h, 
p = resistivity, h = thickness) of this thin resistive plume 
can hardly be resolved into p and h. This normally re- 
sults in smearing with blurred boundaries and larger 
volume relative to the input model. These common 
ERT limitations are minimized here by applying the 
constrained inversion approach taking use of any avail- 
able subsurface data. Uncertainties in mapping struc- 
tures and quantifying the resistivity amplitudes are 
relatively low reflecting the high reliability of the recon- 
structed results. 

Notably we could quantify reliably the gas saturations 
indirectly from the density and resistivity models result- 
ing from the inversion by applying common petrophys- 
ical equations. The saturation results deduced from 
ERT technique fit well their corresponding values de- 
rived from elastic FWI. Both show reasonable absolute 
average differences (<20%) relative to the background. 
However, such results should be cautiously treated, 
where their validity and uncertainty should be studied 
in real field data. 

The use of synthetic data contaminated with random 
error may reflect the real world of data. Obviously, add- 
ing random noise (typically 3%) to the synthetic ERT 
data increases the rms-error values by a factor of 5-9 
but slightly decreases the mapping resolution. Our re- 
sults here are in accordance with that obtained by al 
Hagrey (2012a) for ERT applications in CCS modelling. 
Using modelling codes (as applied here) and adding a 
random noise in an ascending order (1, 2, and 5% levels) 
to the synthetic data sets generally increases the rms- 
errors by a factor of 2 to 9 but slightly decreases the 
mapping capability of ERT technique. Ramirez et al. 
(2005) obtained similar results and concluded that the 
effect of the random error in ERT is insignificant for 
anomalies of a large size and magnitude. 

Obviously the elastic FWI and ERT modelling using 
2.5D codes has been conducted along a 2D section of 
the geological 3D model applied in the modelling simu- 
lation. This 2D model simplification is fully justified by 
the evidence that this 2D section cuts the main (storage) 
structure (gentle anticline within almost horizontal layer- 
ing) along its main strike. 

In conclusion results reveal the capability of our ap- 
plied integrative geophysical approach to resolve the 



CAES targets and to quantify intrinsic property changes 
of the injected gas saturation in the reservoir. Con- 
strained inversion models of elastic FWI and ERT are 
even able to recover well the gradual desaturation with 
depth. The accurately mapped spatial (seismic and elec- 
tric) parameters are applied in their respective petro- 
physical equation to yield precise quantifications of gas 
saturations from each technique independently. Both 
resulting saturation models are in accordance with each 
other and with the input (true) saturation model. A 
joint elastic FWI and ERT inversion has a high poten- 
tial to improve the applicability of the approach (e.g. 
Karaoulis et al. 2012). In a numerical study of a moving 
gas front within a reservoir, the joint inversion of seis- 
mic and electric time-lapse data sets reduces the pres- 
ence of artifacts, and can retrieve the shape and 
estimate parameter better than their individual (uncon- 
strained) inversions. For time-lapse data ERT uses per- 
manently installed borehole electrodes, whereas seismic 
data needs to be repeated such that the source and re- 
ceiver positions could not be the same and therefore 
more uncertainties can occur. 

Moreover, the applied 3D gravity technique shows 
high sensitivity to the mass deficit resulting from the 
storage of the gas phase. The vertical gravity component 
can resolve saturations and saturation changes down to 
±3% assuming that the data is corrected for temporal 
fluctuation effects of the groundwater table. 
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